library(kableExtra)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(readxl)
library(GGally)
library(corrplot)
library(factoextra)
library(ggcorrplot)
library(ggrepel) #para etiquetas
M=matrix(c(3,1,1,1,3,1,1,1,5),nrow=3,byrow=TRUE)
M
## [,1] [,2] [,3]
## [1,] 3 1 1
## [2,] 1 3 1
## [3,] 1 1 5
# Calcular los autovalores y autovectores
eigen_values <- eigen(M)
# Mostrar los resultados
autovalores = eigen_values$values # Autovalores
autovalores
## [1] 6 3 2
autovectores = eigen_values$vectors # Autovectores
autovectores
## [,1] [,2] [,3]
## [1,] -0.4082483 -0.5773503 7.071068e-01
## [2,] -0.4082483 -0.5773503 -7.071068e-01
## [3,] -0.8164966 0.5773503 -1.347436e-16
\(Y_1=-0.4082483 X_1 - 0.4082483 X_2 -0.8164966 X_3\)
\(Y_2= -0.5773503 X_1-0.5773503 X_2 +0.5773503 X_3\)
\(Y_3= 0.7071068 X_1 -0.7071068 X_2 - X_3-1.347436e-16\)
Cada autovalor de la matriz de covarianza en un análisis de componentes principales representa la varianza (o la cantidad de variabilidad) capturada por su correspondiente componente principal. En el contexto de PCA, la suma total de todos los autovalores es igual a la varianza total de los datos originales. La proporción de la variabilidad total de los datos explicada por un componente principal específico se determina dividiendo el autovalor de ese componente principal por la suma de todos los autovalores.
En términos matemáticos, si \(\lambda_1, \lambda_2, \ldots, \lambda_n\) son los autovalores de la matriz de covarianza, entonces la proporción de la variabilidad explicada por el i-ésimo componente principal se calcula como:
\[ \frac{\lambda_i}{\sum_{j=1}^{n} \lambda_j} \]
data.frame(componentes = c('c1','c2','c3'),
Variabilidad = autovalores/sum(autovalores)) %>%
mutate(Variabilidad_acc = cumsum(Variabilidad))
## componentes Variabilidad Variabilidad_acc
## 1 c1 0.5454545 0.5454545
## 2 c2 0.2727273 0.8181818
## 3 c3 0.1818182 1.0000000
data.frame(componentes = c(1,2,3),
Variabilidad = autovalores/sum(autovalores)) %>%
mutate(loading_acc = cumsum(Variabilidad)) %>%
ggplot(aes(x= componentes, y = Variabilidad))+
geom_point()+geom_line()+theme_bw()
autovectores[,1] #loadings de la Componente 1
## [1] -0.4082483 -0.4082483 -0.8164966
data.frame(loading = autovectores[,1], orden = seq(1,3,1))
## loading orden
## 1 -0.4082483 1
## 2 -0.4082483 2
## 3 -0.8164966 3
data.frame(loading = autovectores[,1], orden = seq(1,3,1)) %>%
ggplot(aes(x=orden, y = loading, fill=as.factor(orden)))+
geom_bar(stat = "identity", color='black') +
theme_minimal() +
labs(title = "Loadings de la CP1", x = "Orden", y = "Loading")+
theme_bw()+ scale_fill_brewer(palette = "PuOr")+
theme(legend.position = "none")
X = c(2,2,1)
X
## [1] 2 2 1
Y1= -0.4082483 * 2 -0.4082483 * 2 -0.8164966 * 1
Y2= -0.5773503 * 2 -0.5773503 * 2 + 0.5773503 * 1
scores = as.vector((autovectores[,1:2] * X) %>% colSums())
scores
## [1] -2.449490 -1.732051
for (i in seq(1:2) ){
print(paste0('Y', i, '=', scores[i]))
}
## [1] "Y1=-2.44948974278318"
## [1] "Y2=-1.73205080756888"
Considerando los datos de la base chalets.xls, se pide:
chalets <- read_excel("chalets.xls")
colnames(chalets) = c('promotoraId', 'duracion', 'precio', 'superficie')
chalets = chalets %>% mutate( promotoraId = as.factor(promotoraId))
chalets
## # A tibble: 10 × 4
## promotoraId duracion precio superficie
## <fct> <dbl> <dbl> <dbl>
## 1 1 8.7 0.3 3.1
## 2 2 14.3 0.9 7.4
## 3 3 18.9 1.8 9
## 4 4 19 0.8 9.4
## 5 5 20.5 0.9 8.3
## 6 6 14.7 1.1 7.6
## 7 7 18.8 2.5 12.6
## 8 8 37.3 2.7 18.1
## 9 9 12.6 1.3 5.9
## 10 10 25.7 3.4 15.9
chalets %>%
pivot_longer(cols=2:4, names_to='Variable', values_to='Valor') %>%
ggplot(aes(y=Variable, x = Valor, fill=Variable))+
geom_boxplot()+
scale_fill_brewer(palette = "PuOr")+
geom_jitter(alpha=0.5, shape=21, color='black')+
geom_point(stat = 'summary', fun.data = mean_se, shape=12)+
scale_color_brewer(palette = "PuOr")+
theme_bw()+
facet_wrap(~Variable, ncol=1, scales = "free")
Variables Escaladas:
chalets %>%
mutate_if(is.numeric, scale) %>%
pivot_longer(cols=2:4, names_to='Variable', values_to='Valor') %>%
ggplot(aes(y=Variable, x = Valor, fill=Variable))+
geom_boxplot()+
scale_fill_brewer(palette = "PuOr")+
geom_jitter(alpha=0.5, shape=21, color='black')+
#geom_point(stat = 'summary', fun.x = 'mean')+
scale_color_brewer(palette = "PuOr")+
theme_bw()
chalets[,2:4] %>% summary()
## duracion precio superficie
## Min. : 8.70 Min. :0.300 Min. : 3.10
## 1st Qu.:14.40 1st Qu.:0.900 1st Qu.: 7.45
## Median :18.85 Median :1.200 Median : 8.65
## Mean :19.05 Mean :1.570 Mean : 9.73
## 3rd Qu.:20.12 3rd Qu.:2.325 3rd Qu.:11.80
## Max. :37.30 Max. :3.400 Max. :18.10
Funcion para Valores Atipicos
select_column <- function(dataframe, column_name) {
if (column_name %in% names(dataframe)) {
a = dataframe[[column_name]]
q3 = quantile(a, 0.75)[[1]]
q1 = quantile(a, 0.25)[[1]]
riq = q3-q1
li_salvaje = (q1 - 1.5 * riq)
ls_salvaje = (q3 + 1.5 * riq )
li_severo = (q1 - 3 * riq)
ls_severo = (q3 + 3 * riq )
df = data.frame(Lim_outliers = c('Li_Salvaje','Ls_Salvaje', 'Li_Severo','Ls_Severo'),
valores = c(li_salvaje, ls_salvaje, li_severo, ls_severo ))
tolerance = 1e-6
df2 = dataframe %>% mutate(Categoria = case_when(a < li_salvaje | a > ls_salvaje + tolerance ~ "Salvaje",
a < li_severo | a > ls_severo + tolerance ~ "Severo",
TRUE ~ "Normal")) %>%
filter(Categoria != 'Normal') %>%
#select(column_name, Categoria)
select(all_of(column_name), Categoria)
result = list(valores_criticos = df,
df_filtrado = df2)
return(result)
} else {
stop("El nombre de la columna proporcionado no se encuentra en el dataframe.")
}
}
select_column(chalets, 'duracion')$valores_criticos
## Lim_outliers valores
## 1 Li_Salvaje 5.8125
## 2 Ls_Salvaje 28.7125
## 3 Li_Severo -2.7750
## 4 Ls_Severo 37.3000
select_column(chalets, 'duracion')$df_filtrado
## # A tibble: 1 × 2
## duracion Categoria
## <dbl> <chr>
## 1 37.3 Salvaje
select_column(chalets, 'precio')$valores_criticos
## Lim_outliers valores
## 1 Li_Salvaje -1.2375
## 2 Ls_Salvaje 4.4625
## 3 Li_Severo -3.3750
## 4 Ls_Severo 6.6000
select_column(chalets, 'precio')$df_filtrado
## # A tibble: 0 × 2
## # ℹ 2 variables: precio <dbl>, Categoria <chr>
select_column(chalets, 'superficie')$valores_criticos
## Lim_outliers valores
## 1 Li_Salvaje 0.925
## 2 Ls_Salvaje 18.325
## 3 Li_Severo -5.600
## 4 Ls_Severo 24.850
select_column(chalets, 'superficie')$df_filtrado
## # A tibble: 0 × 2
## # ℹ 2 variables: superficie <dbl>, Categoria <chr>
chalets[,2:4] %>% ggpairs()
Se aprecia correlación entre todos los pares de variables
chalets[,2:4] %>%
summarise_all(mean)
## # A tibble: 1 × 3
## duracion precio superficie
## <dbl> <dbl> <dbl>
## 1 19.0 1.57 9.73
(chalets[,2:4] %>%
summarise_all(mean)) %>%
t() %>%
as.vector()
## [1] 19.05 1.57 9.73
cov(chalets[,2:4])
## duracion precio superficie
## duracion 63.29833 5.7450000 33.863889
## precio 5.74500 0.9934444 4.053222
## superficie 33.86389 4.0532222 20.849000
Tamaño del Problema = Suma de Varianzas de cada Variable.
diag(cov(chalets[,2:4]))
## duracion precio superficie
## 63.2983333 0.9934444 20.8490000
sum(diag(cov(chalets[,2:4])))
## [1] 85.14078
cor(chalets[,2:4])
## duracion precio superficie
## duracion 1.0000000 0.7244728 0.9321763
## precio 0.7244728 1.0000000 0.8906068
## superficie 0.9321763 0.8906068 1.0000000
ggcorrplot(cor(chalets[,2:4]), hc.order = TRUE, outline.col = "white", lab = 'TRUE')
Coeficiente de Determinación.
cor(chalets[,2:4])**2
## duracion precio superficie
## duracion 1.0000000 0.5248609 0.8689527
## precio 0.5248609 1.0000000 0.7931805
## superficie 0.8689527 0.7931805 1.0000000
Si, se trata de variables con correlaciones moderadas a elevadas. Ademas, si elevamos al cuadrado los coef. de correlación obtenemos una aproximación al coeficiente de determinación R2 el cual tambien da valores elevados.
acp = prcomp(chalets[,2:4])
acp
## Standard deviations (1, .., p=3):
## [1] 9.0905247 1.5484743 0.3245999
##
## Rotation (n x k) = (3 x 3):
## PC1 PC2 PC3
## duracion 0.87136808 0.4797725 -0.1026449
## precio 0.08530129 -0.3541685 -0.9312832
## superficie 0.48315770 -0.8027347 0.3495364
El ACP contiene 5 elementos:
sdev: Este elemento contiene las desviaciones estándar de las componentes principales. La longitud de este vector es igual al número de componentes principales. En PCA, las desviaciones estándar se relacionan con la cantidad de varianza que cada componente principal captura del conjunto de datos. Al elevar al cuadrado estos valores, se obtienen los autovalores de la matriz de covarianza o correlación utilizada en el PCA.
rotation: También conocido como la matriz de carga, este elemento contiene los autovectores de la matriz de covarianza/correlación de los datos. Las columnas de esta matriz son los autovectores que corresponden a las componentes principales y que se utilizan para transformar los datos originales al espacio de las componentes principales.
center: Al realizar PCA, es común centrar los datos restando la media de cada variable. Este elemento contiene las medias de las variables originales que se han restado durante el proceso de centrado. Si el argumento center = FALSE se utilizó en prcomp(), este componente no estará presente.
scale:Si los datos fueron escalados (normalmente dividiendo por la desviación estándar de cada variable) antes de aplicar PCA, este elemento contendrá los factores de escala utilizados. El escalado se realiza para que todas las variables contribuyan por igual al análisis, independientemente de sus unidades de medida. Si el argumento scale. = FALSE se utilizó en prcomp(), este componente no estará presente.
x: Contiene los datos originales proyectados en el espacio de las componentes principales, a veces denominados como puntuaciones de las componentes principales (scores). Cada columna es una componente principal y cada fila corresponde a una observación de los datos originales.
acp %>% str()
## List of 5
## $ sdev : num [1:3] 9.091 1.548 0.325
## $ rotation: num [1:3, 1:3] 0.8714 0.0853 0.4832 0.4798 -0.3542 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "duracion" "precio" "superficie"
## .. ..$ : chr [1:3] "PC1" "PC2" "PC3"
## $ center : Named num [1:3] 19.05 1.57 9.73
## ..- attr(*, "names")= chr [1:3] "duracion" "precio" "superficie"
## $ scale : logi FALSE
## $ x : num [1:10, 1:3] -12.33 -5.322 -0.464 -0.269 0.515 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:3] "PC1" "PC2" "PC3"
## - attr(*, "class")= chr "prcomp"
print('sdev')
## [1] "sdev"
acp$sdev
## [1] 9.0905247 1.5484743 0.3245999
print('rotation')
## [1] "rotation"
acp$rotation
## PC1 PC2 PC3
## duracion 0.87136808 0.4797725 -0.1026449
## precio 0.08530129 -0.3541685 -0.9312832
## superficie 0.48315770 -0.8027347 0.3495364
print('center')
## [1] "center"
acp$center
## duracion precio superficie
## 19.05 1.57 9.73
print('scale')
## [1] "scale"
acp$scale
## [1] FALSE
print('x')
## [1] "x"
acp$x
## PC1 PC2 PC3
## [1,] -12.3303278 0.8062792 -7.232126e-02
## [2,] -5.3219077 -0.1712548 2.971035e-01
## [3,] -0.4637910 0.4325717 -4.539599e-01
## [4,] -0.2686924 0.5136236 6.068733e-01
## [5,] 0.5154163 2.0808737 -2.471241e-02
## [6,] -4.8596686 -0.2107264 1.396961e-01
## [7,] 1.2481508 -2.7531684 1.627372e-01
## [8,] 20.0428878 1.6367490 -8.489517e-07
## [9,] -7.4938494 0.0755665 -4.252179e-01
## [10,] 8.9317821 -2.4105140 -2.301978e-01
Para corroborar el analisis, la suma de los valores propios deberia ser igual a la suma de las varianzas. ACP devuelve los desvios de cada Componente por lo cual es necesario elevarlo al cuadrado.
acp$sdev**2 %>% sum()
## [1] 85.14078
acp$rotation
## PC1 PC2 PC3
## duracion 0.87136808 0.4797725 -0.1026449
## precio 0.08530129 -0.3541685 -0.9312832
## superficie 0.48315770 -0.8027347 0.3495364
acp$rotation[,1] %>%
t() %>%
as.data.frame() %>%
pivot_longer(cols = 1:3, values_to = 'LoadingCP1', names_to = 'variable') %>%
ggplot(aes(x=variable, y=LoadingCP1, fill=variable))+
geom_bar(stat='identity', color='black')+
scale_fill_brewer(palette = "PuOr")+
theme_bw()
fviz_eig(acp, addlabels = T, ylim=c(0,120))
acp %>% summary()
## Importance of components:
## PC1 PC2 PC3
## Standard deviation 9.0905 1.54847 0.32460
## Proportion of Variance 0.9706 0.02816 0.00124
## Cumulative Proportion 0.9706 0.99876 1.00000
acp$rotation
## PC1 PC2 PC3
## duracion 0.87136808 0.4797725 -0.1026449
## precio 0.08530129 -0.3541685 -0.9312832
## superficie 0.48315770 -0.8027347 0.3495364
acp$x
## PC1 PC2 PC3
## [1,] -12.3303278 0.8062792 -7.232126e-02
## [2,] -5.3219077 -0.1712548 2.971035e-01
## [3,] -0.4637910 0.4325717 -4.539599e-01
## [4,] -0.2686924 0.5136236 6.068733e-01
## [5,] 0.5154163 2.0808737 -2.471241e-02
## [6,] -4.8596686 -0.2107264 1.396961e-01
## [7,] 1.2481508 -2.7531684 1.627372e-01
## [8,] 20.0428878 1.6367490 -8.489517e-07
## [9,] -7.4938494 0.0755665 -4.252179e-01
## [10,] 8.9317821 -2.4105140 -2.301978e-01
La componente principal 1 (PC1) logra explicar el 97.06% de la variabilidad total de los datos. Al ser todos los loadings positivos se trata de una componente de tamañano.
fviz_pca_biplot(acp, ylim=c(-4,4))
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
Para ordenar las promotoras segun la primer componente seria:
data.frame(id=seq(1:10), cp1 =acp$x[,1]) %>%
arrange(cp1)
## id cp1
## 1 1 -12.3303278
## 2 9 -7.4938494
## 3 2 -5.3219077
## 4 6 -4.8596686
## 5 3 -0.4637910
## 6 4 -0.2686924
## 7 5 0.5154163
## 8 7 1.2481508
## 9 10 8.9317821
## 10 8 20.0428878
data.frame(id=seq(1:10), cp1 =acp$x[,1]) %>%
arrange(cp1) %>% ggplot(aes(x= cp1, y=0, label=id))+
geom_point()+
geom_text_repel(aes(label=id))+
geom_vline(xintercept = c(-4,-8,5), linetype = 2,
color = 2)+
theme_bw()
Es posible agrupar las promotoras según lo manifestado por la CP1. La promotoras 10 y 8 se encuentras separadas del resto. Esto quiere decir que en terminos de superficie y duracion estas promotoras se destacan con clientes de propiedades grandes e hipotecas largas. De manera contraria la promotora 1 se encuentra alejada del comportamiento promedio en terminos negativos, con propiedades chicas y duraciones cortas. Por ultimo, las promotoras 3,4,5,7 son promotoras promedio.
Valores medios de las Variables
acp$center
## duracion precio superficie
## 19.05 1.57 9.73
Para analizar el comportamiento bivariado se realizó un scaterplot entre las variables menos correlacionadas (superficie y duracion). Se puede observar que solo las promotoras 10 y 8 estan por encima de la media en ambas variables. Esto ratifica lo observado en la componente principal1.
chalets %>% ggplot(aes(x=superficie, y=duracion))+
geom_point()+
geom_text_repel(aes(label=promotoraId))+
theme_bw()+
geom_vline(xintercept = 9.73, linetype = 2,
color = 2)+
geom_hline(yintercept = 19.05, linetype = 2,
color = 2)
Dado el siguiente conjunto de datos: \(X = \begin{pmatrix} 3 & 6 \\ 5 & 6 \\ 10 & 12 \\ \end{pmatrix}\)
X = matrix(c(3, 6, 5, 6, 10, 12), nrow = 3, byrow = TRUE)
X
## [,1] [,2]
## [1,] 3 6
## [2,] 5 6
## [3,] 10 12
cov(X)
## [,1] [,2]
## [1,] 13 12
## [2,] 12 12
Tamaño Total del Problema
diag(cov(X)) %>% sum()
## [1] 25
eigen(cov(X))
## eigen() decomposition
## $values
## [1] 24.5104121 0.4895879
##
## $vectors
## [,1] [,2]
## [1,] -0.7216753 0.6922317
## [2,] -0.6922317 -0.7216753
autovalores = eigen(cov(X))$values
autovalores
## [1] 24.5104121 0.4895879
Tamaño Total del Problema
autovalores %>% sum()
## [1] 25
autovectores = eigen(cov(X))$vectors
autovectores
## [,1] [,2]
## [1,] -0.7216753 0.6922317
## [2,] -0.6922317 -0.7216753
data.frame(componentes = c('c1','c2'),
Variabilidad = autovalores/sum(autovalores)) %>%
mutate(Variabilidad_acc = cumsum(Variabilidad)*100)
## componentes Variabilidad Variabilidad_acc
## 1 c1 0.98041649 98.04165
## 2 c2 0.01958351 100.00000
data.frame(componentes = c('c1','c2'),
Variabilidad = autovalores/sum(autovalores)) %>%
mutate(Variabilidad_acc = cumsum(Variabilidad)) %>%
ggplot(aes(x= componentes, y = Variabilidad))+
geom_text_repel(aes(label=round(Variabilidad,2)*100))+
geom_bar(stat='identity', aes(fill=componentes), color='black')+
scale_fill_brewer(palette = 'PuOr')+
theme_bw()
data.frame(X1 = X[,1], X2 = X[,2], id=c(1,2,3)) %>%
ggplot(aes(x=X1, y = X2))+
geom_point()+
geom_text_repel(aes(label=id))+
theme_bw()+
ylim(c(0,12))+
xlim(c(0,13))
Los Scores se pueden realizar de la siguiente manera, multiplicando las
matrices.
scores = X %*% autovectores
scores
## [,1] [,2]
## [1,] -6.318416 -2.2533568
## [2,] -7.761767 -0.8688934
## [3,] -15.523534 -1.7377869
Al aplicar la transformación se pueden diferenciar los puntos 1 y 2 en el sentido de las y, ya que en el original no podria hacerse ya que ambos tenian en valor 6 en Y.
data.frame(cp1 = scores[,1], cp2 = scores[,2], id=c(1,2,3))%>%
ggplot(aes(x=cp1, y = cp2))+
geom_text_repel(aes(label=id))+
geom_point()+
theme_bw()
Xsc = scale(X)
Xsc
## [,1] [,2]
## [1,] -0.8320503 -0.5773503
## [2,] -0.2773501 -0.5773503
## [3,] 1.1094004 1.1547005
## attr(,"scaled:center")
## [1] 6 8
## attr(,"scaled:scale")
## [1] 3.605551 3.464102
cov(Xsc)
## [,1] [,2]
## [1,] 1.0000000 0.9607689
## [2,] 0.9607689 1.0000000
Tamaño Total del Problema: da diferente.
diag(cov(Xsc)) %>% sum()
## [1] 2
eigen(cov(Xsc))
## eigen() decomposition
## $values
## [1] 1.96076892 0.03923108
##
## $vectors
## [,1] [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068 0.7071068
autovalores_sc = eigen(cov(Xsc))$values
autovalores_sc
## [1] 1.96076892 0.03923108
Tamaño Total del Problema
autovalores_sc %>% sum()
## [1] 2
autovectores_sc = eigen(cov(Xsc))$vectors
autovectores_sc
## [,1] [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068 0.7071068
data.frame(componentes = c('c1','c2'),
Variabilidad = autovalores_sc/sum(autovalores_sc)) %>%
mutate(Variabilidad_acc = cumsum(Variabilidad)*100)
## componentes Variabilidad Variabilidad_acc
## 1 c1 0.98038446 98.03845
## 2 c2 0.01961554 100.00000
Los Scores se pueden realizar de la siguiente manera, multiplicando las matrices.
scores_sc = Xsc %*% autovectores_sc
scores_sc
## [,1] [,2]
## [1,] -0.9965967 0.18010011
## [2,] -0.6043644 -0.21213216
## [3,] 1.6009611 0.03203204
data.frame(cp1 = scores_sc[,1], cp2 = scores_sc[,2], id=c(1,2,3))%>%
ggplot(aes(x=cp1, y = cp2))+
geom_text_repel(aes(label=id))+
geom_point()+
theme_bw()
# Tenemos dos autovectores v1 y v2
v1 <- autovectores[,1]
v2 <- autovectores[,2]
vectores_df <- data.frame(
x = c(0, v1[1], 0, v2[1]),
y = c(0, v1[2], 0, v2[2]),
vector = factor(rep(c("Autovector 1", "Autovector 2"), each = 2))
)
vectores_df
## x y vector
## 1 0.0000000 0.0000000 Autovector 1
## 2 -0.7216753 -0.6922317 Autovector 1
## 3 0.0000000 0.0000000 Autovector 2
## 4 0.6922317 -0.7216753 Autovector 2
# Crear el gráfico con ggplot2
ggplot() +
geom_segment(data = vectores_df, aes(x = 0, y = 0, xend = x, yend = y, colour = vector),
arrow = arrow(type = "closed", length = unit(0.2, "inches"))) +
scale_color_manual(values = c("Autovector 1" = "blue", "Autovector 2" = "red")) +
theme_minimal() +
labs(title = "Autovectores en el Espacio Bidimensional Sin Escalar", x = "X", y = "Y") +
coord_fixed()
v1 <- autovectores_sc[,1]
v2 <- autovectores_sc[,2]
vectores_df <- data.frame(
x = c(0, v1[1], 0, v2[1]),
y = c(0, v1[2], 0, v2[2]),
vector = factor(rep(c("Autovector 1", "Autovector 2"), each = 2))
)
vectores_df
## x y vector
## 1 0.0000000 0.0000000 Autovector 1
## 2 0.7071068 0.7071068 Autovector 1
## 3 0.0000000 0.0000000 Autovector 2
## 4 -0.7071068 0.7071068 Autovector 2
# Crear el gráfico con ggplot2
ggplot() +
geom_segment(data = vectores_df, aes(x = 0, y = 0, xend = x, yend = y, colour = vector),
arrow = arrow(type = "closed", length = unit(0.2, "inches"))) +
scale_color_manual(values = c("Autovector 1" = "blue", "Autovector 2" = "red")) +
theme_minimal() +
labs(title = "Autovectores en el Espacio Bidimensional Datos Escalados", x = "X", y = "Y") +
coord_fixed()
La razón por la que los autovectores y autovalores resultan diferentes cuando se calculan con datos escalados y sin escalar radica en cómo la escala de los datos afecta a la matriz de covarianza o correlación que se utiliza para obtenerlos.
Datos sin escalar Cuando realizas un análisis de componentes principales (PCA) sin escalar los datos, estás trabajando con la varianza y covarianza “naturales” de tus datos. Si las variables tienen diferentes unidades o rangos de varianza, las variables con mayor varianza tendrán un peso más significativo en el resultado del PCA. Esto significa que las componentes principales resultantes pueden estar sesgadas hacia estas variables.
Datos escalados Al escalar los datos, normalmente se centran (restando la media) y se estandarizan (dividiendo por la desviación estándar). Este proceso lleva a todas las variables a una misma escala y elimina las unidades, lo que permite que cada variable contribuya equitativamente al análisis del PCA, independientemente de su varianza original o unidades de medida. En otras palabras, el escalado asegura que la varianza de una variable no influya más que otra simplemente debido a la magnitud de sus valores.
Efecto en Autovalores y Autovectores Los autovalores reflejan la varianza capturada por cada componente principal. Si las variables no están escaladas, las componentes principales reflejarán más la varianza de las variables con rangos más grandes. Cuando las variables están escaladas, los autovalores reflejarán una contribución más equitativa.
Los autovectores, que indican la dirección de las componentes principales, también cambian. En datos no escalados, las direcciones pueden estar dominadas por variables con mayor varianza. En datos escalados, cada variable tiene igual oportunidad de influir en la dirección de las componentes principales.
Conclusión La elección de escalar o no los datos antes de realizar el PCA depende del contexto y de los objetivos del análisis. Si las variables están en diferentes unidades o tienen diferentes rangos de varianza y deseas que todas tengan la misma importancia en el análisis, debes escalar. Si algunas variables son intrínsecamente más variables y eso es importante para el análisis, podrías optar por no escalar.
Sea S la matriz de varianzas y covarianzas poblacionales: \(X = \begin{pmatrix} 3 & 1 & 1 \\ 1 & 4 & 0 \\ 1 & 0 & 2 \\ \end{pmatrix}\)
S = matrix(c(3, 1, 1, 1, 4, 0, 1, 0, 2), nrow = 3, byrow = TRUE)
S
## [,1] [,2] [,3]
## [1,] 3 1 1
## [2,] 1 4 0
## [3,] 1 0 2
Crear una matriz de correlacion a partir de una matriz de covarianzas
# Calcular las desviaciones estándar de las variables originales
std_devs <- sqrt(diag(S))
# Crear una matriz de desviaciones estándar para la división
std_dev_matrix <- outer(std_devs, std_devs, "*")
# Calcular la matriz de correlación
cor_matrix <- S / std_dev_matrix
# Mostrar la matriz de correlación
print(cor_matrix)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.2886751 0.4082483
## [2,] 0.2886751 1.0000000 0.0000000
## [3,] 0.4082483 0.0000000 1.0000000
Sí, es posible convertir una matriz de covarianza en una matriz de correlación utilizando operaciones matemáticas estándar y la matriz identidad, junto con la matriz de desviaciones estándar de las variables. La idea es normalizar la matriz de covarianza de tal manera que las varianzas (los elementos de la diagonal) se conviertan en 1, lo cual es característico de una matriz de correlación.
En el contexto de matrices, la diagonal de la matriz de covarianza contiene las varianzas de las variables, y para convertirla en una matriz de correlación, debemos dividir cada elemento de la matriz de covarianza por el producto de las desviaciones estándar asociadas a sus respectivas filas y columnas.
Si D es una matriz diagonal donde cada elemento diagonal D[i, i] es la desviación estándar de la variable i, entonces la matriz de correlación R se puede obtener de la siguiente manera:
\(R=D^{−1}⋅C⋅D^{−1}\) donde C es la matriz de covarianza y \(D^{-1}\) es la inversa de la matriz D (que es simplemente 1 dividido por cada elemento en la diagonal de D, ya que D es una matriz diagonal).
# La matriz diagonal D con las desviaciones estándar
D <- diag(sqrt(diag(S)))
# Invertir D para obtener D^-1
D_inv <- solve(D)
# Calcular la matriz de correlación
cor_matrix <- D_inv %*% S %*% D_inv
# Mostrar la matriz de correlación
print(cor_matrix)
## [,1] [,2] [,3]
## [1,] 1.0000000 0.2886751 0.4082483
## [2,] 0.2886751 1.0000000 0.0000000
## [3,] 0.4082483 0.0000000 1.0000000
autovalores = eigen(S)$values
autovectores = eigen(S)$vectors
autovalores
## [1] 4.732051 3.000000 1.267949
autovectores
## [,1] [,2] [,3]
## [1,] -0.5773503 0.5773503 0.5773503
## [2,] -0.7886751 -0.5773503 -0.2113249
## [3,] -0.2113249 0.5773503 -0.7886751
autovalores %>% sum()
## [1] 9
diag(S) %>% sum()
## [1] 9
autovectores[,2]
## [1] 0.5773503 -0.5773503 0.5773503
Varianza explicada:
autovalores/sum(autovalores)*100
## [1] 52.57834 33.33333 14.08832
autovalores[2]/sum(autovalores)*100
## [1] 33.33333
Para el analisis es necesario recordar de que se tratan las variables.
X1 (Econometría): El peso positivo (0.5744) de la puntuación media obtenida en las asignaturas de econometría sugiere que hay una asociación directa entre las puntuaciones altas en econometría y los valores altos en la segunda componente principal. Es decir, a medida que las puntuaciones en econometría aumentan, también lo hace el valor de la segunda componente principal.
X2 (Derecho): El peso negativo (-0.5744) indica que las puntuaciones altas en derecho están asociadas inversamente con la segunda componente principal. Un valor más alto en derecho contribuiría a un valor más bajo en esta componente principal. Esto podría sugerir que los estudiantes que obtienen altas puntuaciones en derecho tienden a obtener puntuaciones más bajas en econometría y asignaturas libres, o viceversa.
X3 (Asignaturas Libres): Al igual que con X1, el peso positivo (0.5744) significa que las puntuaciones más altas en asignaturas libres se asocian con valores más altos en la segunda componente principal.
La igualdad de los valores absolutos de los pesos indica que la contribución de cada asignatura a la variabilidad capturada por esta componente principal es igual en magnitud. Además, el hecho de que los pesos de X1 y X3 sean positivos y el de X2 sea negativo sugiere que esta componente podría estar capturando un contraste en el desempeño académico: estudiantes que se desempeñan de manera similar en econometría y asignaturas libres pero de manera opuesta en derecho.
En resumen, la segunda componente principal podría estar identificando un patrón de desempeño estudiantil donde existe una dicotomía entre el rendimiento en econometría y asignaturas libres frente al rendimiento en derecho. Estudiantes con altos valores en esta componente principal tenderían a tener puntuaciones similares en econometría y asignaturas libres, pero bajas en derecho, y aquellos con bajos valores en la componente principal mostrarían el patrón opuesto.
Dirección y Ponderación: Las tres variables originales contribuyen de manera igual en magnitud a esta componente principal, ya que los valores absolutos de las entradas del autovector son iguales.
Signo: El signo de las entradas del autovector indica la dirección de la relación. En este caso, la segunda entrada es negativa, mientras que la primera y la tercera son positivas. Esto significa que la segunda variable está inversamente relacionada con la primera y la tercera en esta componente principal.
Interpretación de la Componente: La segunda componente principal parece capturar una variabilidad en los datos donde la primera y la tercera variables varían juntas en la dirección opuesta a la segunda variable. Cabe destacar que la primer variable esta correlacionada positivamente con las variables 2 y 3. Mientras que la segunda variable tiene 0 correlacion con la tercera.
Varianza Explicada: La componente explica el 33% de la variacion de los datos y la igualdad de los valores en el autovector sugiere que esta componente está equidistantemente orientada con respecto a las tres variables.
Relevancia de la Componente: Asumiendo que esta es la segunda componente principal en una PCA, podría representar la segunda mayor fuente de variabilidad en el conjunto de datos después de la primera componente principal.
En resumen, la segunda componente principal parece representar una dimensión de los datos donde hay un equilibrio de variabilidad entre las tres variables observadas, con la segunda variable moviéndose en la dirección opuesta a las otras dos. Para una interpretación más detallada y significativa, necesitaríamos entender qué representan estas variables en el contexto de los datos analizados.
Econometría (X1): Dado que la carga asociada con la econometría es positiva, una puntuación baja en la segunda componente principal sugiere que el estudiante probablemente tiene una puntuación baja en econometría, o al menos más baja en comparación con la puntuación promedio de sus compañeros.
Derecho (X2): La carga negativa asociada con las asignaturas de derecho implica que una puntuación baja en la segunda componente principal indica una puntuación alta en derecho. En otras palabras, el estudiante puede estar destacando en derecho en comparación con sus compañeros.
Asignaturas Libres (X3): Al igual que con la econometría, una carga positiva para las asignaturas libres significa que una puntuación baja en esta componente principal sugiere que el estudiante también tiene un rendimiento bajo en asignaturas libres.
Dado el patrón de cargas en el autovector de la segunda componente principal, una puntuación baja para un estudiante indicaría un rendimiento académico que contrasta con las asignaturas de economía y libres en comparación con el derecho. Es decir, este estudiante podría tener un perfil académico fuerte en derecho pero débil en econometría y asignaturas libres, lo que refleja una inclinación hacia el área legal sobre las otras áreas.
Esto destaca la importancia del contexto en la interpretación de las componentes principales. La segunda componente principal parece capturar una dimensión de variabilidad en el rendimiento académico que distingue entre el enfoque en el derecho y otras áreas de estudio. Por lo tanto, una puntuación baja en esta componente no es necesariamente indicativa de un rendimiento académico general pobre, sino más bien de una distribución específica de fortalezas y debilidades académicas.
cumsum(autovalores/sum(autovalores)*100)
## [1] 52.57834 85.91168 100.00000
data.frame(componentes = c(1,2,3),
Variabilidad = autovalores/sum(autovalores)) %>%
mutate(loading_acc = cumsum(Variabilidad))
## componentes Variabilidad loading_acc
## 1 1 0.5257834 0.5257834
## 2 2 0.3333333 0.8591168
## 3 3 0.1408832 1.0000000
Para obtener al menos el 80% harian falta las componentes 1 y 2.
data.frame(componentes = c(1,2,3),
Variabilidad = autovalores/sum(autovalores)) %>%
mutate(loading_acc = cumsum(Variabilidad)) %>%
ggplot(aes(x= componentes, y = Variabilidad))+
geom_point()+geom_line()+theme_bw()
La tabla gorriones.xls contiene datos de 49 aves, 21 de los cuales sobrevi- vieron a una tormenta. Se pide:
gorriones <- read_excel("gorriones.xlsx")
gorriones = gorriones %>% mutate(sobrevida = case_when(sobrevida == 1 ~ "Si",
sobrevida == -1 ~"No"),
sobrevida = as.factor(sobrevida))
gorriones
## # A tibble: 49 × 7
## pajaro largototal extension cabeza humero esternon sobrevida
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 1 156 245 31.6 18.5 20.5 Si
## 2 2 154 240 30.4 17.9 19.6 Si
## 3 3 153 240 31 18.4 20.6 Si
## 4 4 153 236 30.9 17.7 20.2 Si
## 5 5 155 243 31.5 18.6 20.3 Si
## 6 6 163 247 32 19 20.9 Si
## 7 7 157 238 30.9 18.4 20.2 Si
## 8 8 155 239 32.8 18.6 21.2 Si
## 9 9 164 248 32.7 19.1 21.1 Si
## 10 10 158 238 31 18.8 22 Si
## # ℹ 39 more rows
gorriones_sc= gorriones %>% select(-c(pajaro, sobrevida)) %>%
mutate_if(is.numeric, ~scale(.))
gorriones_sc
## # A tibble: 49 × 5
## largototal[,1] extension[,1] cabeza[,1] humero[,1] esternon[,1]
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -0.542 0.725 0.177 0.0938 -0.329
## 2 -1.09 -0.262 -1.33 -1.01 -1.24
## 3 -1.36 -0.262 -0.578 -0.0900 -0.229
## 4 -1.36 -1.05 -0.704 -1.38 -0.632
## 5 -0.815 0.330 0.0514 0.278 -0.531
## 6 1.37 1.12 0.680 1.01 0.0741
## 7 -0.268 -0.656 -0.704 -0.0900 -0.632
## 8 -0.815 -0.459 1.69 0.278 0.377
## 9 1.65 1.32 1.56 1.20 0.276
## 10 0.00558 -0.656 -0.578 0.645 1.18
## # ℹ 39 more rows
Boxplots Estandarizados
gorriones[,2:6] %>%
mutate_if(is.numeric, scale) %>%
pivot_longer(cols=1:5, names_to='Variable', values_to='Valor') %>%
ggplot(aes(x=Variable, y = Valor, fill=Variable))+
geom_boxplot()+
scale_fill_brewer(palette = "PuOr")+
geom_jitter(alpha=0.5, shape=21, color='black')+
#geom_point(stat = 'summary', fun.x = 'mean')+
scale_color_brewer(palette = "PuOr")+
theme_bw()
Boxplots Sin Estandarizar
gorriones %>%
pivot_longer(cols=2:6, names_to='Variable', values_to='Valor') %>%
ggplot(aes(x=Variable, y = Valor, fill=Variable))+
geom_boxplot()+
scale_fill_brewer(palette = "PuOr")+
geom_jitter(alpha=0.5, shape=21, color='black')+
geom_point(stat = 'summary', fun.data = mean_se, shape=12)+
scale_color_brewer(palette = "PuOr")+
theme_bw()+
facet_wrap(~Variable, ncol=5, scales = "free")
cov(gorriones_sc) %>% round(3)
## largototal extension cabeza humero esternon
## largototal 1.000 0.735 0.662 0.680 0.605
## extension 0.735 1.000 0.674 0.769 0.529
## cabeza 0.662 0.674 1.000 0.756 0.526
## humero 0.680 0.769 0.756 1.000 0.607
## esternon 0.605 0.529 0.526 0.607 1.000
cor(gorriones[,2:6]) %>% round(3)
## largototal extension cabeza humero esternon
## largototal 1.000 0.735 0.662 0.680 0.605
## extension 0.735 1.000 0.674 0.769 0.529
## cabeza 0.662 0.674 1.000 0.756 0.526
## humero 0.680 0.769 0.756 1.000 0.607
## esternon 0.605 0.529 0.526 0.607 1.000
cov(gorriones_sc) %>% round(3) == cor(gorriones[,2:6]) %>% round(3)
## largototal extension cabeza humero esternon
## largototal TRUE TRUE TRUE TRUE TRUE
## extension TRUE TRUE TRUE TRUE TRUE
## cabeza TRUE TRUE TRUE TRUE TRUE
## humero TRUE TRUE TRUE TRUE TRUE
## esternon TRUE TRUE TRUE TRUE TRUE
all.equal(cov(gorriones_sc) %>% round(3) , cor(gorriones[,2:6]) %>% round(3))
## [1] TRUE
En el contexto del Análisis de Componentes Principales (PCA), un autovalor asociado a una componente principal indica la cantidad de varianza que dicha componente principal captura de los datos originales. En términos más técnicos, un autovalor refleja la proporción de la varianza total que es explicada por su componente principal correspondiente.
Magnitud de Varianza: Un autovalor grande indica que esa componente principal captura una gran parte de la varianza en los datos. Por el contrario, un autovalor pequeño sugiere que la componente principal correspondiente capta solo una pequeña cantidad de la varianza.
Orden de Importancia: Las componentes principales se ordenan típicamente por la magnitud de sus autovalores, de mayor a menor. La primera componente principal es la dirección en el espacio de características en la que los datos varían más, seguida de la segunda componente principal, y así sucesivamente.
Elección de Componentes: En muchas aplicaciones de PCA, solo las componentes principales con los autovalores más altos se retienen para reducir la dimensionalidad de los datos. Las componentes con autovalores bajos a menudo se descartan porque capturan menos información (varianza).
Total de Varianza Explicada: La suma de todos los autovalores de un PCA es igual a la suma de las varianzas de todas las variables originales. Esto se debe a que PCA es una transformación ortogonal que conserva la varianza total de los datos.
Contribución Porcentual: La contribución porcentual de cada componente principal a la varianza total se calcula dividiendo el autovalor de la componente por la suma total de los autovalores.
Decisión sobre la Reducción de Dimensiones: A menudo se utiliza un “criterio de corte” como el gráfico de sedimentación (scree plot) o el criterio de Kaiser (mantener solo componentes con autovalores mayores que 1) para decidir cuántas componentes principales retener. También se puede utilizar el porcentaje acumulado de varianza explicada como criterio.
gorriones[,2:6] %>% ggpairs()
acp_gorriones = prcomp(gorriones_sc)
acp_gorriones
## Standard deviations (1, .., p=5):
## [1] 1.9048003 0.7281813 0.6061272 0.5396719 0.4276117
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## largototal 0.4549833 -0.0178440 -0.6172046 0.55597238 0.3203513
## extension 0.4602276 0.3055232 -0.3950893 -0.49111858 -0.5454844
## cabeza 0.4478919 0.3129793 0.5969224 0.49154765 -0.3217162
## humero 0.4727851 0.1642198 0.2941304 -0.44718825 0.6804529
## esternon 0.3962927 -0.8839757 0.1418940 -0.09000453 -0.1824950
summary(acp_gorriones)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.9048 0.7282 0.60613 0.53967 0.42761
## Proportion of Variance 0.7257 0.1061 0.07348 0.05825 0.03657
## Cumulative Proportion 0.7257 0.8317 0.90518 0.96343 1.00000
Hacen falta las dos primeras componentes para representar al menos el 80% de la variabilidad
fviz_eig(acp_gorriones, addlabels = T, ylim=c(0,100))
Totas las CP.
cp_gorriones = acp_gorriones$rotation
cp_gorriones
## PC1 PC2 PC3 PC4 PC5
## largototal 0.4549833 -0.0178440 -0.6172046 0.55597238 0.3203513
## extension 0.4602276 0.3055232 -0.3950893 -0.49111858 -0.5454844
## cabeza 0.4478919 0.3129793 0.5969224 0.49154765 -0.3217162
## humero 0.4727851 0.1642198 0.2941304 -0.44718825 0.6804529
## esternon 0.3962927 -0.8839757 0.1418940 -0.09000453 -0.1824950
La primer componente:
cp_gorriones[,1]
## largototal extension cabeza humero esternon
## 0.4549833 0.4602276 0.4478919 0.4727851 0.3962927
\(CP1 = largototal * 0.4549833 + extension * 0.4602276 + cabeza * 0.4478919 + humero * 0.4727851 + esternon * 0.3962927\)
gorriones %>% filter(pajaro==11)
## # A tibble: 1 × 7
## pajaro largototal extension cabeza humero esternon sobrevida
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 11 158 240 31.3 18.6 22 Si
p11 = gorriones %>%
filter(pajaro==11) %>%
select(-c(pajaro, sobrevida))%>%
t() %>%
as.vector()
p11
## [1] 158.0 240.0 31.3 18.6 22.0
CP1:
cp_gorriones[,1]
## largototal extension cabeza humero esternon
## 0.4549833 0.4602276 0.4478919 0.4727851 0.3962927
cp_gorriones[,1] %*% p11
## [,1]
## [1,] 213.8732
p11[1] * 0.4549833 + p11[2] * 0.4602276 + p11[3] * 0.4478919 + p11[4] * 0.4727851 + p11[5] * 0.3962927
## [1] 213.8732
fviz_pca_biplot(acp_gorriones, geom = c("point", "text"),
col.ind = "cos2", # Color por calidad de la representación
col.var = "contrib", # Color por contribución a la PCA
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Evitar solapamiento de texto
)
fviz_pca_biplot(acp_gorriones, axes = c(1, 3), geom = c("point", "text"),
col.ind = "cos2", # Color por calidad de la representación
col.var = "contrib", # Color por contribución a la PCA
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Evitar solapamiento de texto
)
fviz_pca_biplot(acp_gorriones, axes = c(2, 3), geom = c("point", "text"),
col.ind = "cos2", # Color por calidad de la representación
col.var = "contrib", # Color por contribución a la PCA
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Evitar solapamiento de texto
)
cbind(acp_gorriones$x, gorriones$sobrevida) %>%
data.frame() %>%
select(PC1, PC2, V6) %>%
mutate(sobrevida = case_when(V6 == 2 ~ "Si", V6 == 1 ~"No") ) %>%
ggplot(aes(x=PC1, y=PC2, color=sobrevida))+
geom_point()+
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Línea horizontal en y = 0
geom_vline(xintercept = 0, linetype = "dashed", color = "black") + # Línea vertical en x = 0
theme_bw()